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We present a new, simple renormalization group method of investigating groundstate properties of 
interacting bosonic systems. Our method reduces the number of particles in a system, which makes 
numerical calculations possible for large systems. It is conceptually simple and easy to implement, 
and allows to investigate the properties unavailable through mean field approximations, such as 
one- and two-particle reduced density matrices of the groundstate. As an example, we model a 
weakly interacting ID Bose gas in a harmonic trap. Compared to the mean-field Gross-Pitaevskii 
approximation, our method provides a more accurate description of the groundstate one-particle 
density matrix. We have also obtained the Hall-Post lower bounds for the groundstate energy of 
the gas. All results have been obtained by the straightforward numerical diagonalization of the 
Hamiltonian matrix. 

PACS numbers: 03.65.-w, 05.10.Cc, 67.40.Db 
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I. INTRODUCTION 



The numerical investigation of the groundstate prop- 
erties of a multiparticle interacting bosonic system is 
a much harder task than in the case of a single par- 
ticle system. The naive approach consists in choosing 
a large enough finite Hilbert space basis and the nu- 
merical diagonalization of the resulting Hamiltonian ma- 
trix. However, the necessary basis size grows exponen- 
tially with the number of particles, which makes this 
simple method inadequate for the treatment of large sys- 
tems. To avoid this problem, many approximations have 
been invented, such as the Gross-Pitaevskii (GP) mean 
field approach [l], the Density Matrix Renormaliza- 
tion Group (DMRG) method 0, i] or, in the case of 
strong interactions, the Thomas-Fermi approximation [l[ 
and the Tonks-Girardeau model 0, [||. For fermionic 
systems, the Exact Diagonalization- Ab Initio (EDABI) 
method has been implemented @, H, One should 
also note the exceptional case of a full analytical solu- 
tion in one dimension by Lieb and Liniger [1 01 ] - From 
this solution, two- and three-pair correlation functions 
of an interacting one-dimensional (ID) Bose gas have 
been derived [H E3, El- In this paper we present a 
new approach, which is in the flavour of Renormalization 
Group methods but is conceptually simple and easy to 
implement. Our method amends the problem of unman- 
ageable basis size by reducing the number of particles in 
the system and renormalizing the Hamiltonian. We ap- 
proximate one- and two-particle properties of the large 



system with the same properties of the smaller system. 
In contrast to mean field methods, our approach allows to 
calculate such quantities as one- or two-particle reduced 
density matrices (1-RDM and 2-RDM) of the ground- 
state. 

In Section HH we describe how our method works. Sec- 
tion IIIII contains an example application of the method 
to the problem of one-dimensional interacting Bose gas 
in a harmonic trap. The results are summarized in Sec- 
tion hyi 



II. HAMILTONIAN RENORMALIZATION AND 
THE APPROXIMATION OF GROUNDSTATE 
PROPERTIES 

We present our approximation in the case of a system 
with two-body interactions. It can be easily generalized 
to the general case of n-body interactions. 

Consider a ID Hamiltonian describing a system of N 
scalar (zero spin) bosons, 



N 

E 

k=\ 



h 2 d 2 

2m dx\ 



+ Vi(aJk) 



N k-1 
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E V l( X k' X k>) , 
k=l k' = l 

(I) 

where N is the number of particles, m is the particle 
mass, Vi(x) is the external one-particle potential and 
V2(x,x') — V2{x',x) is the two-particle interaction po- 
tential. 

Investigating such a system, we are often interested 
only in one- or two-particle properties of the groundstate. 
One way to calculate them is to obtain an approximation 
of the 1-RDM or 2-RDM of the groundstate. 

We approximate the system by replacing it with a 



2 



much smaller system containing N' <§C N scalar bosons. 
The smaller system is governed by a renormalized version 
of the original Hamiltonian Hn-, that is 



H N . = — 
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iV 7 
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2m dxi 
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N' fe-1 



T—[Y1 E V 2( x k,x k ') 



fe=l fc': 



(2) 



We calculate the properties of interest from the ground- 
state of ffjv' j thus avoiding the insurmountable problem 
of diagonalizing the Hamiltonian of the large system, H^. 
The results for increasing values of N' will converge to the 
values of the corresponding properties of the ./V-particle 
system. 

We will now justify our procedure. Let \1/ and p be a 
state of the large system and its iV'-RDM, respectively. 
It is easy to show that their mean energies, measured by 
the respective Hamiltonians, are equal, 



Tr [H N ,p \ = (* H 



IN 



(3) 



Hence, when the mean energy of ^ becomes lower, mov- 
ing closer to the mean energy of the groundstate \&o of 
Hn, the mean energy of p, also becomes lower and moves 
closer to the mean energy of the (pure state) density ma- 
trix p' of the groundstate ty' of Hn'- Because of the 
variational principle, the density matrix p' is an approx- 
imation of the reduced density matrix (RDM) po of the 
groundstate ^o- The one- or two-particle properties of 
po (i.e. of ^o), like the probability density, are approx- 
imated by the same properties of p' (i.e. of *S? ). Since 
N' <C N, it is much easier to calculate numerically the 
groundstate than the groundstate ^o, and to investi- 
gate the one- or two- particle properties of "Jo by inves- 
tigating the same properties of 

The main source of error in our method is the fact 
that the variational search for the groundstate con- 
verges to the ./V'-particle groundstate 'J'q, not to the 
RDM of the TV-particle groundstate, pq. This is be- 
cause for bosons, not every iV'-particle density matrix 
is an RDM of an iV-particle state. A better strat- 
egy would be to perform the variational search not in 
the whole iV'-particle Hilbert space but in the smaller 
space of iV'-particle density matrices which are RDMs 
of iV-particle states. However, the problem of identi- 
fyin g th is space, the so-called iV-representability prob- 
lem [H HE El, remains unsolved. Therefore, we have to 
perform our calculations for a sequence of N' . The energy 
of the iV'-particle groundstate is a lower bound of the 
energy of the investigated iV-particle groundstate [l7l ]. 
When N' increases, the groundstate energy increases 
and approaches the groundstate energy of the ./V-particle 
groundstate. Due to the variational principle, this means 
that the ./V'-particle groundstates approximate the N- 
particle groundstate increasingly well, in the sense that 
the one- and two-particle properties calculated from these 



groundstates converge to the corresponding properties of 
the iV-particle groundstate. 

In the general case of n-particle inter- 
actions, the renormalization goes as fol- 
lows: an rt-particle interaction potential term 

J2k ± =i J2k 2 =i ' ' ' J2kl=i ' ' ' J2k n =i Vn(%i, ■ ■ ■ , x n ), 
symmetrical with respect to permutations of coordi- 
nates Xk, is replaced by the term 



Efc=iEfe=i'"Eit°=i v n (xi, 

is true also in this general case. 



jyny-i) (iv-n+i) 

N'(N'-l) (iV'-n+l) 

,x n ). Equation © 



III. A SIMPLE EXAMPLE 

A. Investigated system 

In our example, we consider a system of N = 100 scalar 
bosons with a dimcnsionless Hamiltonian 



fe=i 



Xk') 



k=l k' = l 



where 5{x 



is the Dirac 5 function. This interaction 



potential is often used to describe cold bosons forming 
a Bosc-Einstcin condensate (BEC) in a trap, when only 
the s-wave scattering occurs Our example concerns 
positive A, which lead to the repulsive interaction. 

We have approximated numerically the 1-RDM and 
2-RDM of the groundstate for two values of interaction 
strength A. The procedure begins with the calculation of 
a finite matrix of the renormalized Hamiltonian Hj^i in a 
finite basis composed of the noninteracting Hamiltonian 
(A = 0) eigenstates, permanents 19] of N' one-particle 
Hermite functions Ti.k, 



H k {x) 



1 



=H k (x)exp \^- — 



(4) 



where H^ix) is fc-th Hermite polynomial. The basis 
contains all eigenstates with (nonrcnormalized) energies 
lower than a cutoff energy L + N' /2, i.e. permanents 

of functions 7ik n , n = 1, . . . , N' such that J^Li < L. 
Then, the groundstate is calculated with the help of an it- 
erative Lanczos-type numerical procedure (l8| . From the 
groundstate we obtain the 1-RDM and the 2-RDM, with 
trace normalized to unity. Using them, we can calculate 
any one- or two-particle property of the groundstate. For 
given N', the basis cutoff L is chosen to be large enough 
so that calculated properties do not change upon further 
increase of L. 

In the case of 1-RDM, p%, we compare the diagonal 
part, pi(x,x), with the probability density calculated by 
minimizing numerically the GP energy functional [l( of 



3 



our system, 



©[*] 



N 
~2 



dx 2 



+ x 2 ^(x) 



+ (JV-l)A|*(a;)r 



(5) 



dx . 



The minimization is performed by expanding the wave- 
function ^(x) in the finite basis of the first twenty Her- 
mite functions ([3]), inserting the expansion into ([5]) and 
minimizing numerically the resulting functional of the 
expansion coefficients. 

We present numerical results for two values of A, 10~ 2 
and 5 • 10 -2 . Both values are in the regime of weak inter- 
actions, where the minimization of the GP energy func- 
tional provides a good approximation of the groundstate. 



B. Groundstate energy 

First, we provide the data on the groundstate energy 
Eq. In the noninteracting case A = 0, Eq is precisely 
known and equals 50. Table U lists three different ap- 
proximations of Eq for two nonzero values of A. In the 
second and third column of Table Q] two different upper 
bounds for Eq are listed: the one obtained from the GP 
functional, ©gp, and the variational bound, ©Gauss, cal- 
culated as a minimal mean value of Hjq in a state a 
product of N Gaussian one-particle wavefunctions with 
a common variational parameter a, 



*<r(a;i 



N 1 



fe=i 



: exp 



4<7 2 



i.e. ©Gauss = min (TeR (* .|i/jv| , I , <T)- Relatively small 
differences between ©gp and -©Gauss indicate that the 
groundstates are close to Gaussian. As expected, the 
GP approximation provides a better estimation of the 
groundstate energy than the Gaussian Ansatz. Our 
method provides an estimation of the true groundstate 
energy Eq as the so-called Hall-Post [ijj lower bound, 
Enp. Values of ©hp are listed in the fourth column of 
Table HJ They were calculated by NLSQ (nonlinear least 



A 


Eqp 


^Gauss 


Enp 


10 -2 
5 • 1(T 2 


68.8 
130.9 


68.9 
132.3 


68.2 
121.4 



our calculations, we have used results for N' = 8, so as 
to make ©hp as high as possible) as function of L to a 
power law 

E(L) = ©hp + BL C , C < 

and taking the L — > oo limit, obtaining ©hp as the an- 
swer. It has been necessary to follow this procedure, 
since raw numerical results vary with L, even for L large 
enough so that the 1-RDM does not change. The rela- 
tive asymptotic standard error of the fitting procedure is 
below 0.01% for both values of A. 

The bounds on Eq (Enp and ©gp) presented above are 
quite close. The relative uncertainty with which ©o is de- 
termined by them is below 1% for A = 10~ 2 and around 
8% for A = 5 • 10 -2 . (We have calculated the relative un- 
certainty as the ratio of the difference between the upper 
and the lower bound to the lower bound.) The fact that 
it is small supports the applicability of our approxima- 
tion, as it means that the true groundstate energy is also 
close to the obtained lower bound. On the other hand, 
if the reduced density matrix of the iV'-particle ground- 
state is to be a good approximation of the reduced 
density matrix of the true groundstate, the mean energy 
of \Pq — i.e. the Hall-Post lower bound — must be close to 
the true groundstate energy of the system. Our results 
fulfill this condition. 



C. Density matrices 

For A = 10 -2 , we obtain identical one-particle prob- 
ability densities from our method and from the GP ap- 
proximation, as shown on Fig. [TJ The accuracy of our 



N'=5, L=50 
GP 
nonint. 




Table I: Bounds for the groundstate energy Eq. The first two, 
Egp and -EGauss, are the upper bounds calculated from the 
GP approximation and the variational method with a Gaus- 
sian wavefunction, respectively. The last one, -Ehp, is the 
Hall-Post lower bound, which is provided by our method as 
an estimate of the true groundstate energy. 

squares) fitting of the groundstate energy for fixed N' (in 



Figure 1: Comparison of the probability densities obtained 
from our method (for N' — 5 and L = 50) and from the GP 
approximation, for A = 10 -2 . The curves overlap perfectly, 
indicating convergence. The probability density of the nonin- 
teracting (A = 0) groundstate is shown too. One can clearly 
see the difference between the interacting system and the non- 
interacting one, with the repulsive interaction pushing away 
the bosons in the trap. 



4 



approximation is confirmed by Fig. [2j which shows that 
different values of N' and L yield identical probability 
densities. A magnified section of this plot is shown in 




-3-2-10123 



x 

Figure 2: Probability densities for A = 1CP 2 calculated for 
increasing N' (L = 56, 50, 60, 40, 40 respectively). The curves 
overlap perfectly, indicating convergence. 

Fig. [3J We will use the convergence with increasing N' 
as a benchmark of the accuracy of our method, treating 
our numerical results as correct if they stabilize quickly. 
For each N' , we take the results for L large enough so 
that they do not change upon further increase of L. A 
similar convergence occurs for the antidiagonal part of 
the f-RDM, pi(x,-x). 

The GP approximation, however, cannot provide us 
with the knowledge about the nondiagonal parts of the 
f-RDM. The merit of our method is that we can calcu- 
late px(x,y) for any (x,y). For A = 10~ 2 , we obtain 
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Figure 3: Probability densities for A = 10~ 2 calculated for 
increasing N' (L = 56,50,60,40,40 respectively), shown in 
the range x £ [—0.25, 0.25]. The plot has been maginified in 
order to show the details of the convergence. 



numerically 



pi(x,y) ~ pi [vs 2 +y 2 J , (6) 

which is clearly shown by the contour plot of pi(x,y) in 
Fig.Hl 
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Figure 4: Contour plot of pi(x,y) for A = 10 2 , calculated 
for N' = 5 and L = 50. Isolines display a radial symmetry. 

The convergence of the diagonal part of 1-RDM, 
pi{x, x) and of the diagonal part of 2-RDM, p2{x, x, x, x) 
(the probability of finding both particles in the same po- 
sition a;), is compared in Fig. [51 ft is clear that the con- 
vergence, with increasing N', of the second function is 
much slower. The consequence of this is that using only 
such simple diagonalization techiques as we did, which 
limit us to N' < 10, we cannot estimate the 2-RDM, 
even for A as small as 10 ~ 2 . 



0.25 




Figure 5: On the left plot, the probability density pi(x, x) for 
A = 10~ 2 is shown for increasing N' (L — 56, 50, 60, respec- 
tively). The curves converge quickly. On the right plot, the 
diagonal part of the 2-RDM, p<z{x, x, x, x) is plotted (for the 
same A) for increasing TV' (L — 56, 50, 60, 40, 40, respectively). 
Even for N' equal 7 or 8, the curves do not converge. 
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For A = 5- 10 -2 , we obtain the convergence of the prob- 
ability densities as easily, as for A = 10~ 2 (see Fig. [S]), 
although it is slightly slower (not visible on the plot). 
A magnified section of this plot is shown in Fig. [7J A 




Figure 6: Convergence of the probability densities pi(x, x) for 
A = 5- 10" 2 and increasing values of N' (L = 55, 50, 60, 40, 40, 
respectively). The curves overlap, indicating convergence. 

similar convergence occurs for the antidiagonal part of 
the 1-RDM, pi(x, — x). However, the probability den- 
sity differs slightly from the one obtained from the GP 
functional, as seen in Fig. [H Contrary to the case of 
A = lfr 2 , the contour plot of the 1-RDM for A = 5 • 1(T 2 
is no longer radially symmetric, as can be seen on Fig. [5] 
It differs noticeably from the one (shown in Fig. [TO)) we 
would obtain from the mean field method, using the for- 
mula 




Figure 8: Probability density for A = 5 • 10~ 2 , as calculated 
with our method (TV' = 8 and L = 40) and from the minimiza- 
tion of the GP functional. A slight difference between the two 
curves can be seen in the middle of the plot, indicating that 
interactions are strong enough so that GP approximation's 
results differ from ours. 
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Figure 7: Probability densities for A = 5 • 10~ 2 calculated for 
increasing N' (L — 55,50,60,40,40 respectively), shown in 
the range x G [—0.25,0.25]. The plot has been magnified in 
order to show the details of the convergence. 




Figure 9: Contour plot of pi(x,y) for A = 5- 10 -2 , calculated 
for N' = 8 and L = 40. Isolines do not display the radial 
symmetry, unlike for A = 10~ 2 . 



where ^gp is the real wavefunction which minimizes the 
GP energy functional. This difference is the most striking 
result in this Section, and indicates that our method gives 
more accurate results than mean- field approximations. 

For even higher value of A, 10 -1 , we did not obtain 
fast enough convergence of neither p±(x, x) (see Fig. ITT|) 
nor, especially, p±(x, — x) (see Fig. [T2|) . This prevented 
us from investigating the full 1-RDM for this interaction 
strength. We conclude that A = 1CP 1 is outside the range 
of our approximation in its present form. To reach this 
interaction strength, we would have to calculate the den- 
sity matrices of the iV'-particle groundstate for N' higher 
than those treatable with the simple numerical diagonal- 
ization algorithm used by us. 
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Figure 10: Contour plot ol pi(x, y) calculated for A = 5 ■ 1CF 2 
using GP approximation. It is clearly visible that this plot is 
more symmetrical than the one in Fig. [9] 
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Figure 11: On the left plot, the probability density pi(x,x) 
for A = 10~ 2 is shown for increasing N' (L = 56, 50, 60, re- 
spectively). On the right plot, the the same quantity is shown 
for A = 10" 1 (L = 55,50,60, respectively). The convergence 
for the higher value of A is visibly slower. 



IV. SUMMARY 

We have presented a method of investigating one- 
and two-particle reduced density matrices of the ground- 
state of an interacting systems with a large number of 
bosonic particles. The method approximates it with a 
smaller, renormalized system. It is conceptually sim- 
ple and easy to implement numerically. The results it 
provides would be, for high enough interaction strengths 
(e.g. A = 5 • 10~ 2 in our example), impossible to calcu- 
late using mean-field methods, such as the GP approxi- 



mation. 

We have provided an example application of our 
method to the problem of one-dimensional interacting 
Bose gas in a harmonic trap, obtaining accurate approx- 
imations of the quantity unobtainable from mean field 
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Figure 12: Antidiagonal part of 1-RDM for A = lO^ 1 for in- 
creasing N (L = 55, 50, 60, 39, 40, respectively). Convergence 
is too slow for the results to be meaningful. 



methods, namely the full one-particle density matrix. 
The results are precise and accurately describe the large 
system, which is proven by the fact the the results con- 
verge quickly with increasing N'. The GP approximation 
does not give as accurate picture of the groundstate one- 
particle density matrix as our approach. Additionally, 
the Hall-Post lower bounds for the groundstate energy 
have been calculated. Relatively small difference between 
them and the upper bounds (GP and Gaussian) also sup- 
ports the applicability of our method. 

Even using simple numerical procedures, our method 
gives access to properties which were previously not as ac- 
curately described by mean-field methods. To investigate 
the two-particle density matrix, or to perform simula- 
tions of systems with higher interaction strengths, would 
require the use of more refined approaches to the calcula- 
tion of the groundstate of the renormalized Hamiltonian, 
for example the DMRG method @, H or the variational 
optimization of the basis wavefunctions @, H, || . 
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